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Abstract 

The  empirical  testing  of  a program  often  calls  for  generating  a set  of  random 
numbers  and  then  immediately  sorting  them.  In  this  paper  we  consider  the  problem 
of  accomplishing  that  process  In  a single  step:  generating  a sorted  list  of  random 
numbers  (specifically,  reals  chosen  uniformly  from  [0,1]).  The  method  we  describe 
generates  the  randoms  in  linear  time,  Is  perfectly  random  (if  It  can  call  a perfectly 
random  generator  for  a single  uniform),  and  can  be  described  In  Just  three  lines  of 
Algol  or  Pascal  code.  If  the  numbers  are  not  required  to  be  generated  all  at  once 
(but  are  rather  to  be  used  one-at-a-tlme),  then  the  method  can  be  Implemented  as  a 
subroutine  to  produce  the  ^nexf^riumber  and  requires  only  constant  storage. 
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1.  Introduction 


The  first  step  of  many  computer  algorithms  Is  to  sort  the  Input  data.  When 
testing  these  programs  to  determine  runtimes  empirically,  one  usually  generates  N 
random  numbers  and  then  sorts  them.  The  efficiency  of  sorting  algorithms  Is 
well-known  (see  Knuth  [1973]),  however,  so  It  Is  often  not  necessary  to  test  the 
sorting  procedure  empirically  In  a particular  program.  In  this  application  (as  well  as 
many  others),  It  Is  desirable  to  generate  a sorted  Hst  of  random  numbers  as  quickly 
as  possible.  In  this  paper  we  will  study  the  problem  of  generating  a sorted  list  of  N 
reals  drawn  uniformly  from  [0,1]. 

The  most  obvious  method  for  generating  a sorted  list  of  randoms  Is  to  first 
generate  N randoms  (see  Knuth  [1969,  Chapter  3])  and  then  sort  them.  This 
method  requires  time  proportional  to  N Ig  N In  the  worst  case,  but  this  can  be 
reduced  to  linear  expected  time  if  a "bucket"  sort  Is  used  (see  Knuth  [1973, 
Section  5.2.1  ]).  This  linear  expected  time  algorithm  Is  rather  complicated  to  code, 
and  requires  extra  space  proportional  to  N.  A knowledge  of  elementary  probability 
theory,  however,  allows  one  to  use  more  sophisticated  approaches. 

In  this  paper  we  will  Investigate  a new  method  for  generating  sorted  lists  of 
randoms  that  has  significant  advantages  over  all  previous  approaches.  We  will 
begin  by  discussing  previous  work  In  Section  2.  In  Section  3 we  will  study  some 
Important  probabilistic  lemmas,  and  then  show  In  Section  4 how  these  can  be  used 
to  make  efficient  programs.  A discussion  of  this  approach  Is  offered  In  Section  5. 

2.  Previous  Work 

Before  presenting  our  new  algorithms  for  generating  sorted  lists  of  random 
numbers,  we  will  mention,  for  purposes  of  completeness  and  of  comparison,  the  best 
previously  known  method  for  generating  sorted  lists  of  random  numbers.  Although 
the  method  seems  to  be  well-known  among  statisticians,  the  present  authors  ere 
unable  to  find  a description  of  Its  computational  aspects  In  the  statistical  literature. 
The  algorithm  Is  based  on  the  following  lemma. 
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Lemma  1_: 

If  Xv  X2,  ...  • Xn+1  are  Independent  random  vertablea  with  exponential 
distribution  of  any  fixed  mean,  then  the  values 

Yj-[  Zx,]  / [ Z X,] 

1SKJ  1SKn+1 

for  J ■ 1 , . . . , n are  distributed  as  the  order  statistics  of  a random  sample  of 
size  n from  La[0,1]. 

Proof: 

We  omit  the  proof  of  this  lemma  as  It  Is  well-known  (see,  for  example,  Johnson 
and  Kotz  [1070,  Chapter  18])  and  Is  not  essential  to  the  main  thrust  of  this 
paper.  □ 

An  algorithm  derived  from  Lemma  1 Is  described  by  the  following  pseudo-Pascal 
code.  It  assumes  that  RAND  Is  a function  that  on  each  call  returns  an  Independent 
random  number  from  the  uniform  distribution  on  [0,1];  a random  exponential  Is  then 
achieved  by  negating  the  natural  logarithm  of  RAND.  The  effect  of  the  algorithm  Is 
to  fill  elements  1 ..N  of  the  array  X with  sorted  random  numbers  Independently  drawn 
from  U[0,1  ]. 

Sue  «■  Oi 

for  I *■  1 to  N do 

XII)  * Sum  «•  Sue  - ln(RAND)t 
Sue  * Sue  - ln(RAND)t 
for  I * 1 to  N do 
XU)  «■  XII]/Suet 

Program  1^  Filling  an  array  with  sorted  randoms. 

It  Is  obvious  that  this  method  Is  a very  efficient  way  of  generating  sorted  lists  of 
numbers  chosen  uniformly  on  [0,1].  Its  one  computational  disadvantage,  however,  la 
that  It  Is  Inherently  a two-pass  algorithm— the  first  to  place  the  numbers  Into  the 
array  and  the  second  to  normalize  them.  We  wM  now  turn  our  attention  to  a new, 
single-pass  algorithm. 
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The  correctness  of  the  algorithms  to  be  presented  in  Section  4 rests  on  the 
following  two  lemmas.  Lemma  2 will  allow  us  to  generate,  In  constant  time,  the 
largest  of  n independent  uniformly  distributed  random  numbers.  Lemma  3 shows  that 
once  we  have  generated  the  k largest  of  n Independent  uniform  randoms,  the 
problem  of  generating  the  k+18*  largest  reduces  to  the  problem  of  generating  the 
largest  of  n-k  Independent  uniform  randoms. 

Lemma  2: 

The  probability  distribution  of  the  maximum  of  n Independent  random  numbers 
from  the  distribution  U[0,1]  Is  the  same  as  that  of  the  n^1  root  of  a single 
number  from  U[0,1]. 


Proof: 


Note  that  the  both  distributions  mentioned  range  over  the  Interval  [0,1].  Let 

q € [0,1].  It  suffices  to  show  that  numbers  from  either  distribution  have  equal 
probability  of  being  In  [0,q]. 

If  X Is  drawn  from  U[0,1],  then  P(X^n<q)  * P(X<qn)  « qn.  On  the  other  hand, 
the  largest  of  a set  of  n numbers  In  [0,1]  will  lie  In  [0,q]  Iff  all  n lie  In  [0,q]. 
Since  the  probability  of  a single  number  drawn  from  U[0,1]  will  lie  In  [0,q]  Is  q, 
It  follows  that  the  probability  of  n numbers  drawn  Independently  from  U[0,1  ] all 
being  less  than  q Is  also  qn.  □ 


Lemma  3: 


Let  n and  k be  positive  Integers,  n < k.  Let  yj,  ....  yk  be  elements  of  [0,1] 

such  that  y-|  £ y2  £ . . . * y^.  Then,  for  n random  numbers  X1 Xn  chosen 

Independently  from  U[0,1]  the  distribution,  conditional  on  the  largest  k being 

yj Vk*  of  the  k*18t  largest  Is  the  same  as  the  distribution  of  the  largest 

of  n-k  numbers  uniformly  selected  from  [O.yi.]. 
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Proof1 : 

We  note  first  that  the  probability  of  the  k+18t  largest  of  the  Xj  being  equal  to 
yk  Is  zero,  as  Is  the  probability  of  the  largest  of  n-k  Independent  draws  from 
U[0,yk]  being  equal  to  yk.  It  remains  to  consider  the  case  where  the  k+1 8t 
largest  of  the  Xj  lies  In  [0,yk). 

Consider  the  event  space  of  all  sets  of  n independent  draws  from  U[0,1j.  The 
subspace  containing  all  events  wherein  the  largest  k numbers  drawn  are 
y V • > Yk  maV  be  partitioned  Into  a number  of  equivalence  classes.^  Each 
such  equivalence  class  may  be  obtained  by  assigning  the  y(  to  k of  the  X|, 
picking  all  events  from  the  full  space  which  satisfy  these  assignments,  and 
throwing  away  all  events  In  which  any  of  the  n-k  "unspecified"  X|  happen  to  be 
larger  than  or  equal  to  yk.  Thus,  the  distribution  of  the  smallest  n-k  entries, 
within  each  equivalence  class,  Is  precisely  the  distribution  of  n-k  Independent 
draws  from  U[0,yk).  Since  there  are  finitely  many  equivalence  classes,  It 
follows  that  the  distribution  of  the  n-k  smallest  entries,  within  the  union  of  all 
equivalence  classes  (/.e.,  contingent  only  on  the  k largest  draws  being 
y 1 , • . • , yk  and  on  the  k+18t  being  less  than  yk)  Is  Identical  to  the  distribution 
of  n-k  Independent  draws  from  U[0,yk).  This  completes  consideration  of  the 
case  In  which  the  k*1s1  largest  of  the  X(  Is  In  [0,yk),  so  we  are  done.  O 

4.  Programs 

In  this  section  we  will  see  how  the  basic  probabilistic  facts  discussed  In  the  last 
section  can  be  used  to  make  programs  for  generating  sorted  lists  of  randoms.  In  all 
these  programs  we  will  assume  that  we  have  a subroutine  RAND  that  returns  a 
random  number  drawn  uniformly  from  [0,1].  All  the  programs  that  we  will  describe 
produce  correct  output  In  the  sense  that  If  RAND  satisfies  the  probabilistic  definition 
of  U[0,1],  then  the  output  of  our  program  will  satisfy  the  probabilistic  definitions  of  a 

^ Since  this  paper  is  intended  primarily  for  non-statisticians,  we  have  attempted  to  minimize  statistioal  notation  in 
the  presentation  of  this  lemma,  af  the  expense  of  conciseness.  A more  general  form  of  this  well-known  result  la 
more  formally  presented  as  Theorem  2.7  of  David  [1970]. 

2This  number  (the  number  of  event  classes)  may  range  from  1 to  nl/(n-k)t,  depending  on  the  number  and  pattern 
of  equalities  among  the  y^. 
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sorted  list  of  N such  randoms. 

Lemma  2 of  Section  3 allows  us  to  generate  the  maximum  of  N uniforms  In  [0,1  ] 

by  evaluating  RANDAN,  which  we  will  call  CurMax  (for  reasons  soon  to  become 

obvious).  After  we  have  done  that,  Lemma  3 allows  us  to  solve  the  remainder  of  the 

problem  by  generating  N-1  randoms  uniform  on  [0, CurMax].  We  can  accomplish  this 

by  taking  as  the  maximum  the  value  of  CurMax-RAND1^-1),  and  so  forth.  This 

process  can  be  described  precisely  by  the  following  program  In  pseudo-Pascal, 

which  places  the  random  numbers  Into  the  array  X in  decreasing  order. 

CurMax  * 1.0s 

for  I * N downto  1 do 

XIII  «■  CurMax  «■  CurMax  * RAND*/!  s 

Program  2.  Straightforward  Implementation. 

In  the  above  program  the  variable  CurMax  represents  the  current  maximum  of  the 
range  In  which  I randoms  are  to  be  generated.  (A  program  essentially  equivalent  to 
Program  2 was  described  by  Friedman  [1971]  for  use  In  random  event  generation  In 
a physics  context.  He  did  not,  however,  observe  the  generality  of  his  method.) 

In  Program  2 we  exponentiate  to  a fractional  power.  Since  most  programming 
languages  do  not  directly  support  such  a statement,  this  step  Is  usually  Implemented 
as 

XII]  «■  CurMax  CurMax  * exp(ln(RAND)/I). 

The  multiplication  In  that  statement  might  be  a source  of  numerical  error,  so  It  can 

be  replaced  by  an  addition  as  In  the  following  program  to  fill  X with  sorted  randoms. 

LnCurMax  «*  0,0i 
for  I * N downto  1 do 
begin 

LnCurMax  «•  LnCurMax  ♦ ln(RAN0)/It 
X[I]  * exp(LnCurMax) 
end  i 

Program  3.  Multiplication  replaced  by  addition. 


Note  that  with  perfect  arithmetic  this  procedure  will  produce  exactly  the  same 
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output  as  Program  2 (assuming  the  use  of  the  same  procedure  RAND);  this  program, 
however,  la  numerically  more  robust  than  Its  predecessor. 

In  many  applications  the  variables  are  not  all  needed  at  one  time,  but  rather  can 
be  used  "on  the  fly".  If  this  Is  Indeed  the  case,  then  using  the  N array  words  of  X Is 
very  wasteful  of  storage.  We  would  prefer  to  have  an  algorithm  that  can  generate 
the  "next"  value.  We  will  now  describe  such  an  algorithm  as  two  Pascal 
subroutines.  Procedure  InltSorted  Is  passed  an  Integer  N and  Initializes  the  global 
variables  I (an  Integer)  and  LnCurMax  (a  real);  NextSorted  Is  a parameterless 
function  that  returns  the  next  value  In  decreasing  order  (unless  N values  have 
already  been  returned). 

procedure  InltSorted(Nt  Integer)) 
begin 
I «■  Ns 

LnCurMax  «■  0,0 
ends 

function  NextSorted!  reals 
begin 

If  I <»  0 then  Aborts 
LnCurMax  +>  LnCurMax  ♦ ln(RAND)/Is 
I «■  I-ls 

NextSorted  «■  exp(LnCurMax) 
ends 


Program  4.  On-line  generation  of  sorted  randoms. 

Making  N successive  calls  on  NextSorted  after  executing  InltSorted(N)  will  produce 
exactly  the  same  output  as  executing  either  Program  2 or  Program  3 (although  not 
In  the  array  X).  If  an  N+1st  call  Is  made  on  NextSorted  then  abnormal  termination 
will  be  effected  by  calling  procedure  Abort.  As  this  algorithm  la  stated  It  returns  the 
values  In  decreasing  order;  If  Increasing  order  Is  preferred  then  this  can  be 
accomplished  by  subtracting  the  result  from  one. 
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6.  Discussion 

Programs  3 and  4 of  the  previous  section  have  been  Implemented  as  Pascal 
programs;  these  programs  are  described  by  Bentley  and  Saxe  [1970].  Both 
Implementations  required  approximately  250  microseconds  to  generate  a single 
random  number  when  executed  on  a Digital  Equipment  Corporation  PDP-10  KL 
processor.1  To  compare  these  programs  to  more  straightforward  methods  of  solving 
this  problem  we  wrote  a program  that  generates  an  array  of  N random  uniforms  and 
then  uses  Quicksort  to  sort  the  array.  The  Implementation  of  Program  3 was 
somewhat  slower  than  the  sorting  methods  for  values  up  to  N * 250;  after  that  point 
Program  3 Is  faster.  A significant  advantage  of  our  programs  over  the  naive 
methods,  however,  Is  that  while  the  sorting  algorithm  was  described  by  some  80 
lines  of  Pascal  code,  our  method  requires  only  a dozen  lines.  To  ensure  that  the 
randomness  properties  of  our  algorithms  were  not  adversely  affected  by  roundoff 


errors  or  by  using  a llnear-congruentlal  psuedo-random  number  generator,  we  ran  a 
number  of  statistical  tests  to  determine  the  randomness  of  the  resulting 
numbers--all  tests  were  passed  with  flying  colors. 

Throughout  this  paper  our  programs  have  taken  logarithms  of  real  numbers 
uniformly  distributed  on  [0,1].  Notice  that  this  leads  to  an  undefined  result  If  the 
value  of  the  random  number  Is  zero.  Although  this  does  not  affect  the  theory 
underlying  the  paper  (since  we  only  took  such  logarithms  to  "simulate"  fractional 
exponentiation  or  generate  exponentially  distributed  randoms),  this  will  affect 
programs  Implementing  these  methods.  Such  programs  should  take  logarithms  of 
randoms  uniform  on  (0,1].  Since  most  RAND  subroutines  return  values  uniform  on 
[0,1),  this  can  be  accomplished  by  using  1-RAND  as  the  desired  random  number. 

Although  it  Is  clear  that  the  method  of  Program  3 Is  superior  to  a 

Hhe  Pascal  compiler  used  in  these  tests  does  not  proSice  very  efficient  code*  the  authors  suspect  that  the 
speed  of  the  programs  could  be  sUutantially  increased  by  careful  hand-coding.  This  Is  unnecessary  In  meet 
applications,  however,  since  the  use  of  this  method  is  usually  enough  to  remove  the  process  of  generating  sorted 
randoms  from  the  time  bottleneck  of  the  program. 
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generate-and-sort  solution  In  almost  all  applications,  It  Is  more  difficult  to  compare 
Program  3 with  Program  1 . Program  1 Is  faster  than  Program  3 (Program  1 uses  an 
addition,  a logarithm,  a multiplication,  and  three  array  accesses  for  each  random; 
Program  3 uses  an  additional  exponentiation,  but  only  one  array  access),  but 
Program  3 Is  shorter  to  code.  The  primary  advantage  of  the  method  of  Section  4. 
over  Program  1 Is  that  this  method  can  be  Implemented  on-line;  the  method  of 
Program  1 has  no  on-line  version  corresponding  to  Program  4. 

Although  we  have  described  our  method  for  generating  sorted  lists  of  uniform 
random  numbers,  the  same  method  can  be  extended  to  generate  sorted  numbers 
from  other  distributions.  To  generate  numbers  from  distribution  F for  which  the 
inverse  F"1  Is  known,  it  is  only  necessary  to  apply  the  monotone  function  F"1  to 
each  of  the  outputs  of  Programs  3 or  4,  and  the  resulting  sorted  list  will  satisfy  all 
the  desired  properties. 
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ALGORITHM.  Generating  Sorted  Lists  of  Randoms 

DESCRIPTION 

This  Pascal  program  implements  two  algorithms  described  by  Bentley  and  Saxe 
[1070]  for  generating  sorted  lists  of  random  numbers.  The  theory  underlying  these 
algorithms  can  be  found  In  that  paper. 

REFERENCE 

Bentley,  J.  L.  and  J.  B.  Saxe  [1979].  "Generating  sorted  lists  of  random 
numbers,"  attached. 


ALGORITHM 

(ft  ROUTINES  FOR  GENERATION  OF  SORTED  RANDOM  NUMBERS  ft) 

(ft  The  algorithms  used  in  this  program  are  taken  from  "Generating 
sorted  lists  of  random  numbers",  hereinafter  referred  to  as 
"Bentley  and  Saxe."  The  reader  should  refer  to  that  article  for 
a discussion  and  justification  of  the  algorithms.  The  procedure 
GenSorted  implements  Program  3 of  Bentley  and  Saxe  for  filling 
an  array  with  sorted  random  numbers  uniformly  drawn  from  the 
interval  10,11.  The  procedure  InitSorted  and  the  function 
NextSorted  together  Implement  Program  4 of  Bentley  and  Saxe  for 
generating  sorted  random  numbers  on-line.  The  main  program  ie  a 
test  driver  which  exercises  these  routines.  ft) 


: | 


const 

MaxRands  ■ 100;  (ft  Maximum  number  of  random  numbers 

generated  by  Gensorted  ft) 

TeatSize  - 25»  (ft  Number  of  sorted  randoms  to  generate — 

used  by  test  driver.  (Notes  TeetSize 
must  be  <-  MaxRands)  ft) 
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var 

(ft  Variables  for  on-line  generation  of  sorted  randoms.  These 
variable  names  are  the  same  as  those  used  in  Program  4 of 
Bentley  and  Saxe,  except  that  they  have  been  preceded  by 
"OLG"  (for  On-Line  Generation)  to  lessen  the  probability  of 
name  conflicts  with  other  global  variables  which  may  occur 
in  programs  using  the  on-line  generation  routines.  ft) 

OLGI:  integer?  (*  The  next  random  number  generated  by 

NextSorted  will  be  the  OLGI-th 
smallest,  ft) 

OLGLnCurflax:  real;  (ft  The  natural  logarithm  of  the  previous 

number  generated  by  NextSorted.  Before 
the  first  call  of  a sequence,  LnCurMax 
is  set  to  0,  i.e.,  I n (1 ) . ft) 

(ft  Variables  used  by  driver  ft) 

Jt  integer? 

Y?  RandArray? 

(*  Storage  used  by  underlying  random  number  generator  ft) 

RandHold?  integer? 


(ft  Procedures  for  generation  of  uniform  random  numbers  ft) 

(*  The  built-in  function.  Random,  takes  a single  integer  argument 
and  returns  a pseudo-random  real  number  in  the  range  (0,1).  The 
argument  (here  named  RandHold)  is  a VAR  parameter  used  to  hold 
the  current  random  seed,  and  is  altered  by  each  call  to  Random.  *) 
(ft  Note?  The  function  Random  is  not  a Standard  Pascal  built-in 
function.  At  your  site  the  random  number  function  may  go 
by  a different  name,  or  it  may  even  be  necessary  for  you 
to  urite  your  oun.  (See  CALGO  Algorithms  266  and  294  or 
Section  3.2  of  Knuth’s  The  Art  of  Computer  Programming, 
Volume  2t  Semi -Numerical  Algorithms,  Addleon-Uesley, 
1963.)  Also,  the  method  of  initializing  the  random 
number  generator  may  vary  from  site  to  site.  In  short, 
the  programmer  should  be  prepared  to  rewrite  the  routines 
Rand  and  InitRand  to  conform  to  the  local  runtime 
environment.  ft) 

procedure  InitRand? 
begin 

RandHold  i-  0 
end? 

function  Rand?  real? 
begin 

Rand  i-  1 -Random (RandHold)  (ft  return  a number  in  (0,1).  *) 
end? 
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(*  Routine  to  place  N random  numbers  uniformly  draun  from  [0,11 
into  XI1..N1  in  ascending  order.  *) 

(*  The  algorithm  used  here  is  that  of  Program  3 of  Bentley  and 
Same.  *) 

procedure  GenSortedfvar  X«  RandArrayt  Ni  Integer)) 
var 

It  Integer) 

LnCurflax:  real) 
begin 

LnCurflax  t ■ 0.0) 
for  I !■  N down to  1 do 
begin 

LnCurflax  :»  LnCurflax  + In  (Rand) /I) 

X C 1 1 t-  exp  (LnCurflax) 
end 

end) 


(ft  Routines  to  generate  sorted  randoms  on-line  ft) 

(*  To  generate  N random  numbers  from  (0,1),  sorted  in  descending 
order:  First  call  Ini tSor ted(N) . The  next  N evaluations  of 
NextSorted  will  return  the  random  numbers,  in  descending  order. 

If  InitSort  is  called  again  before  N calls  have  been  made  to 
NextSorted,  the  current  sequence  of  randoms  will  be  lost  and  a 
new  sequence  Mill  begin  with  the  next  call  to  NextSorted.  *) 

(ft  Note:  I f an  ascending  sequence  of  random  numbers  is  desired,  the 
final  assignment  statement  of  Nextsorted  should  be  altered  to 
read  "NextSorted  :-  1 - exp (OLGLnCurflax) " . ft) 

(ft  The  algorithms  used  here  are  from  Program  4 of  Bentley  and  Saxe. 

The  global  variable  names,  I and  LnCurflax,  occurring  in  that 
program  have  here  been  preceded  by  "OLG"  to  guard  against  their 
accidental  use  by  other  pieces  of  code.  ft) 

procedure  Ini tSorted(N:  integer)) 
begin 
OLGI  t-  Nt 
OLGLnCurflax  :■  0.0 
end: 

function  NextSortedt  real) 
begin 

i f OLGI  <-  0 then 
begin 

uriteln(tty,  'Too  many  calls  on  NextSorted,  Aborted')) 

ha  1 1 

end) 

OLGLnCurflax  t - OLGLnCurflax  + I n (Rand) /OLGI  | 

OLGI  t-  OLGI -1 | 

NextSorted  :■  exp  (OLGLnCurflax) 
and) 
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(*  Driver  to  test  both  single-shot  and  on-line  generation  of  sorted 
randoms  *) 


begin  (ft  main  program  ft) 
Ini tRand; 


(*  Test  single-shot  generation  of  sorted  randoms  by  filling 
Ytl.  .TestSize)  with  ascending  sorted  random  numbers  and  dumping 
results  to  output  file.  ft) 

GenSorted(Y,  TestSize); 

ur  i telnfDump  of  array  Y after  execution  of  GenSorted’); 
for  J ;•  1 to  TestSize  do 
uritelnU,  Y[JJ:10:5)j 

(*  Test  on-line  generation  of  sorted  randoms  by  generating  TestSize 
random  numbers  in  descending  order  and  uriting  them  to  output,  ft) 
ur  i te In; 

ur  I tel  n(’  Commencing  test  of  on-line  generation'); 

Ini tSor ted (TestSize) ; 
for  J i*  1 to  TestSize  do 

ur itein(J.  NextSorted;  10; S) 
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